clear all; clc; close

%[(1)200 (2)400 (3)400 (4)400 (5)400 (6)400 (7)600 (8)800 (9)800 (10)1000 (11)1200 (12)2000]

%load EcoliHyperShockEpi.mat
load CompileHyperShockData.mat
expIndex = 3;

[numCells,~]=size(Lengths{expIndex}); % numCells is number of cells


%      R = 0.000555556; % The average unperturbed growth rate
%      deltaT = R.*(t_vec{i});
%      L0 = exp(deltaT);

[~,numPoints]= size(Time{expIndex});

locationOfShock = 21;
     
     %numPoints = locationOfShock; % when plotting only before the jump. 
 
% figure
% hold on
% 
% xlabel('Time (s)')
% ylabel('L/L0')
% title(ShockMags(expIndex))

     
    %% calculating the average length as function of time
     avg_length = zeros(numPoints,1);
     N = zeros(numPoints,1);
     for j=1:numCells
         if (~isnan(Lengths{expIndex}(j,locationOfShock))) && (~(Lengths{expIndex}(j,locationOfShock)==0))
           
           normalized_length = Lengths{expIndex}(j,:)/Lengths{expIndex}(j,locationOfShock);
           
           %plot(Time{expIndex},Lengths{expIndex}(j, :)/Lengths{expIndex}(j, locationOfShock))
           %size(Lengths{expIndex}(j, :))
             
           normalized_length(isnan(normalized_length)) = 10;
           

        %plot(deltaT,normalized_length)
        for k=1:numPoints
            if normalized_length(k) ~= 10
           
         avg_length(k) = (normalized_length(k) +  avg_length(k));
         N(k) = N(k) + 1;
            end
        end
         end
       % plot(deltaT,normalized_length)
     end
     %avg_length
     avg_length = avg_length ./ N;
    % n
     %N
     
%      figure
%      hold on
%      
%      plot(Time{expIndex}, avg_length)
%      xlabel('Time (s)')
%      ylabel('L/L0')
%      title(ShockMags(expIndex))
%      
%      fig2pretty
    

     
        %% calculating the standard deviation corresponding to the mean above
     std_length = zeros(numPoints,1);
     N = zeros(numPoints,1);
     for j=1:numCells
         if (~isnan(Lengths{expIndex}(j,locationOfShock))) && (~(Lengths{expIndex}(j,locationOfShock)==0))
           
           normalized_length = Lengths{expIndex}(j,:)/Lengths{expIndex}(j,locationOfShock);
           
%          normalized_length = cell_length{i}(j,:)/cell_length{i}(j,1);
%          normalized_length = normalized_length./L0 - 1;
           normalized_length(isnan(normalized_length)) = 10;
           

        %plot(deltaT,normalized_length)
        for k=1:numPoints
            if normalized_length(k) ~= 10
           
         std_length(k) = std_length(k) +  (normalized_length(k) -  avg_length(k))^2;
         N(k) = N(k) + 1;
            end
        end
         end
       % plot(deltaT,normalized_length)
     end
     %avg_length
     std_length = sqrt(std_length ./ N);
    % n
     %N
     
%      figure
%      hold on
%      
%      plot(Time{expIndex}, std_length)
%      xlabel('Time (s)')
%      ylabel('std')
%      title(ShockMags(expIndex))
%      
%      fig2pretty
     
     
     
     figure
     hold on
     
     plot(Time{expIndex}, avg_length, 'b-')
     
     plot(Time{expIndex}, avg_length - 2*std_length, 'r-')
     
     plot(Time{expIndex}, avg_length + 2*std_length, 'r-')
     
     xlabel('Time (s)')
     ylabel('L/L0')
     title(ShockMags(expIndex))
     
     hLeg = legend('example');
     set(hLeg,'visible','off')
     
     fig2pretty
    
     
     %% Export arrays
     
     expName = ['LBto' num2str(ShockMags(expIndex))];
     
     writematrix(Time{expIndex},['time_'  expName '.csv']) 
     writematrix(avg_length,['avgLength_' expName '.csv']) 
     writematrix(std_length,['stdLength_' expName '.csv']) 
     
 %%
%      fullSamples = [];
%      for j = 1:numCells
%          
%        if  (~isnan(cell_length{i}(j,locationOfShock))) && (~(cell_length{i}(j,locationOfShock)==0))
%             fullSamples = [fullSamples j];
%        end
%      end      
%      
%      figure 
%      hold on
%      
%      
%      for j=fullSamples
%         
%      end
%    